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1. Introduction 



The production of hadronic jets in high-energy collisions is one of the most striking 
and universal features of particle physics. It is now well understood that jets are as- 
sociated with hard parton emission, followed by parton showering and hadronization. 
However, the observed features of jets depend not only on their intrinsic properties 
but also on the algorithms used to find them. These algorithms have been the object 
of many years of development and refinement - see [1] for a review. 

The jet algorithms used most widely nowadays are of the sequential recombina- 
tion type. Of these, the k t algorithm [2-4] was used extensively at LEP and HERA, 
and to a limited extent at the Tevatron. The good theoretical properties of the k t 
algorithm mean that jet rates and multiplicities can be calculated perturbatively, 
either in fixed order or in an all-orders logarithmic approximation. However, the 
irregular angular shapes of fc^-jets make them less suitable for analyses of hadron 
collider events, where multiple parton interactions give rise to underlying activity 
not associated with the primary hard process. 

A variant of the k t algorithm that allows better control of the underlying event 
is anti-kt [5], which has become the preferred jet finder at the LHC It belongs to 
the general category of inclusive generalized k t jet algorithms, defined in [6]. Many 
of the good theoretical features of k t extend more broadly to the members of this 
category. In particular their leading and next-to- leading double logarithmic jet rates 
can be calculated to all orders, and are in fact the same for all members, including 
kt and anti-fcj. 

In the present paper we calculate these jet rates, along with the average jet 
multiplicities, in a variety of approximations. We do so by finding, in each case, the 
equations that the quark and gluon jet generating functions [7-9] must satisfy. In 
some cases, an explicit solution in terms of special functions can be found, while in 
others we are forced to resort to numerical methods (though we can get some analytic 
understanding of the properties of the solution, as we describe in detail for one case in 
an Appendix). We then compare the results of our theoretical calculations with those 
obtained via Monte Carlo simulations using the SHERPA event generator [10,11]. 
Our theoretical calculations of jet rates are carried out for the case of e + e~ collisions, 
but we also perform calculations and simulations of the average jet multiplicity in 
hadron-hadron collisions, finding that similar behaviour is obtained. 

The outline is as follows. We begin by recalling the definitions of the inclusive 
jet algorithms and the relevant Sudakov form factors. In Section 4, we introduce 
the jet generating functions and derive the integral equations that they satisfy at 
double-leading-logarithmic accuracy (DLA) and next-to-double-leading-logarithmic 
accuracy (NDLA). In Sections 5 and 6, we compute the resulting jet rates and av- 
erage jet multiplicities, both at fixed coupling and also at next-to-leading order in 
the running coupling. In Sections 7 and 8, we compare with results of simulations 
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performed using SHERPA for e + e~ and pp collisions, respectively. In Section 9 we 
consider the implications for the scaling patterns of jet multiplicities. Our conclu- 
sions are summarized in Section 10. Details of the derivation and properties of the 
partial differential equation (PDE) for the average jet multiplicity are relegated to 
appendices. 



2. The inclusive generalized kt jet algorithms 

We consider first the case of multijet production in e + e~ annihilation, for which the 
inclusive algorithms are defined as described in the Fast Jet user manual [6], Sect. 4.5. 
The distance measures are 

d ij =mm{E i »,E j >} {1 _ cqsR) 

= min{£f, Ef , 
d iB = E? , (2.1) 

with p — 1, 0, —1 for the k t [2], Cambridge/Aachen [12,13] and anti-/c 4 [5] algorithms, 
respectively. At any stage of clustering, if a dy is the smallest measure we combine 
objects i and j. If diB is the smallest we call i a jet candidate and remove it from the 
clustering list. We then call jet candidates with energy Ei > Er resolved jets. Thus 
the jet rates, at a given value of the e + e~ centre-of-mass energy E cm , are functions 
of the radius parameter R and the minimum jet energy Er. This is in contrast to 
the exclusive k t (Durham) algorithm, where one effectively sets £r = | and continues 
clustering objects until all d^ are above some fixed value d cut = y C ntE 2 m , so that the 
jet rates are functions of the single parameter y cut . 

In hadron-hadron collisions, the cm. frame of the parton-parton hard scattering 
process is not known and therefore one has to adopt a longitudinally invariant form 
of the algorithms [3,4]. To that end, Eqs. (2.1) are replaced by 

AR 2 

d ij = min fe 2 f^?f}-^ , 
ARl = ( Vi - y 3 f + (0, - , 
diB = Ptf , (2.2) 

where p t i, y% and 4>i are the transverse momentum, rapidity and azimuth of object i, 
respectively, and we define jet candidates with p ti > E R as resolved jets. 

As far as leading logarithms are concerned, the jet rates defined by (2.1) and 
(2.2) will be the same, and therefore in the following Sections we refer mainly to 
e + e~ annihilation. By "leading logarithms" here we always mean leading double and 
next-to-double logarithms, a§ log 2n and a§ log 2 ™ -1 , where the logarithms are those of 
1/R and/or Q/Er, Q being the hard process scale. By taking Er sufficiently large 
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in hadron-hadron collisions, we avoid such leading contributions from initial-state 
showering and the underlying event, so these terms are determined by the timelike 
showering of final-state partons. 



3. Sudakov factors 



The evolution scale for coherent parton showering is £ = 1 — cos 9 with 9 the emission 
angle. To be resolved, an emission must have £ > £r and E > Er. The probability 
for a single resolvable gluon emission from a quark of energy E at scale £ is thus 



(3.1) 



where the running coupling is evaluated at the transverse momentum scale of the 



emission, kf = z E £', 



a s (k 



7i = b \n{z 2 E^'/A 2 ) ^ 
with b = (11CU — 2n/)/12. Defining a s = o:s(E 2 ^)/tt, i.e. in terms of the coupling 
at the hard scale, we have to next-to-double-log accuracy (NDLA) 
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(3.4) 



Then the probability for no resolvable emissions (the quark Sudakov factor) is 

A q (E,Z)=exp[-V q (E,Z)] . (3.5) 
Similarly for a gluon, the probability of a single resolvable gluon, quark or antiquark 



emission is 



v 9 (e,o 

which gives to NDLA 

V g (E,t) = a a ]n(J^ 



Er/E 



dz a ^[P 99 (z) + P q9 (z)} , 



(3.6) 
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(3.7) 
(3.8) 



and the gluon Sudakov factor is 

A a (E,t)= eX p[-r a (E,t)] . 

Note that all this is independent of the value of p, so that all the inclusive generalized 
kt algorithms are equivalent at this level of precision. 
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4. Generating functions 

By definition the generating function for resolved jets from a quark (i = q) or gluon 
(i = g) of energy P at scale £ is [7-9] 

oo 

$ l (u,E,0 = Y. unR n(E,0, (4-1) 

n=0 

where R l n is the corresponding n-jet rate, i.e. the probability of finding n resolved 
jets. Thus the jet rates can be recovered from the generating function by successive 
differentiation at u — 0: 

Rl i (E,0 = ^^Uu,E,0\u=o. (4.2) 

On the other hand the average multiplicity of resolved jets is obtained by differenti- 
ating at u — 1. Writing the average jet multiplicity from a quark or gluon of energy 
E at scale £ as A/i(P, £), we have 

M t (E,0 = J2nK(E,0 = ^Si(ti,£,OI«=i • (4-3) 

n=0 

The generating functions <& qtg must thus satify the boundary condition 

$ t (u, E, = 1 + (u - 1)0(P - E R ) . (4.4) 

The generating function for e + e~ annihilation at cm. energy E cm is that for two 
quarks of energy E cm /2, each filling one hemisphere: 

$ ee = [<I>>,P cm /2,l)] 2 . (4-5) 

4.1 Next to double leading logarithms 

For £ > £ R and E > E R , we have to NDLA 

%(u,E,0=uA q (E,0+ f ^^P 5 ^)*>,P^(^P,0 , 

+P^)[$ (? ( M ,P,0] 2 }. (4.6) 

Using the expressions given earlier to eliminate the Sudakov factors, we may 
write these in the equivalent forms 



$ q (u, P, = u + I ^ f dz°^-P gq (z)$ q (u, P, O zE, O - 1] , 

y^H ? y*v£ 27r 



+P 9S (z) [{$„(«, P, OY - P, O] } • (4-7) 
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The solution for the quark generating function is then easily seen to be 

'U £' Je r /e 



* q {u,E,Q = ue*p\ [ %- f dz°^-P gq (z)[%(u,zE,a-l}) 

Ue K ? Je r IE 27T J 



We can solve for the gluon generating function by iteration, and then substitute in 
this equation to get the complete solution. 

4.2 Leading double logarithms 

In the DLA we keep only the singular parts of P gq and P gg and can drop P qg . For 
brevity we define the logarithms as 

K = ln(E/E R ), A = ]n(£/&)- (4.9) 

Then 

A gi9 (M)=e- a - KA (4.10) 

where a q<g = Cp^as, and 

$ g (u,K, A) = ue~ a,?KA exp \a q J dn' J d\' $ g (u, k', A')| , (4.11) 

$ g (u, K, A) = e- a ° KX ju + a 9 jf d«' jf dA' e a ^ v $ 3 (u, «, A')$ 5 (m, A') 



We can simplify the equation for $ 9 by noting that 

( e a ^ A $ g ) = a g e a ° KX $ g f dn $ g (u, k , A) (4.12) 



so that 



\n(e a ^ x <S> g ) = [ dn' [ d\'$ g (u,Kf,\') + C(u,K) (4.13) 
Jo Jo 

where the boundary condition <& g (u, K, 0) = it for all k > implies C(u,k) = lnw. 
Thus 

$ p (u,k,A) = «e _<vtA exp ja 9 jf da/jf dA' ® g (u, k', A') J • (4.14) 
Comparing with (4.11) we see that in the DLA 

$ g = u(<f> g /uf F/CA . (4.15) 

5. Jet rates 

From the quark and gluon generating functions in the previous section it is relatively 
straightforward to construct the jet rates to next-to-double log accuracy (NDLA). 
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In terms of the logarithmic variables (4.9) the quark and gluon integrated splitting 
kernels (3.4) and (3.7) are simply 



P q (E,£)=a s \C F 



~-C F b al\K [2A + k] 



(5.1) 



V g (E, = a s A [C a k - b ] + -C A b Q al\K [2A + «] 



(5.2) 



where 



a s = (xs(k,\) 



1 



b (2K + A + n) 



(5.3) 



with/x = ]n(££6i/A 3 ). 

The quark generating function is 



$ g (u,E,\) = wexp |c F dA'^ ^a s (« / , A') - |a s e"' -,e ^ [$< 



(u,«', A') - 1] 



(5.4) 



In the fixed-order expansion it is not necessary to include the running coupling for 
the finite terms in the splitting functions, as doing so would affect results beyond 
the NDLA. In these finite term we can set terms proportional to e~ K to after the 
integrations, which is equivalent to allowing the original z integration to range over 
(0, 1), since it is not singular in energy. However, we have to bear in mind that the 
relevant integrals vanish when k — 0. 

Defining 



Tg(n', A', k) = C/ 



T q (K , A , «) = C F 



a s («',A') - ^«se K ' K 



a s («',A') - -a s e K K 



(5.5) 
(5.6) 
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we solve the gluon generating function by iteration to third order in u. 



& g 1 \u,K,\)=uA 9 (K,\) (5.7) 

^ g 2 \ U ,K,X) =uA g (K,\) (l + U [ dW dK? ^T g (K',X',K) A g (K',X') 

(5.8) 



+ ^^l^ 



6 A s (k,A') 



<$>W(u,k,X) = uA 9 (k,X) [ 1 + u I d\' I dK'{r g (K',\',K)A g (K',\') 

-A' 



X 1+U 



[ K +r dK" [ d\" \t 9 (k", A", k)A g {K n , A") 
Jo Jo Jo I 



+ as 6A s (,',A") 6 



+ a s ^^p^-e^(l + 2uC F f* dX" /V 
6 A fl (K,A') V A Jo 

xr,( K ",A", K )A 9 ( K ",A") 



(5.9) 



In the result for & g we have defined R = k(k) when the k integral ranges from 
to k (/«'). Defining 



¥ g n) = 2 I dK' I dX T q (K', X, «) ^ n) (u, A') 
Jo Jo 



(5.10) 



the generating function for resummed rates up to 5 jets at NDLA is then 



$i 5 e ) = u*A%k, A) {l + *W + 1 ) 2 + I 



(5.11) 



Substituting in (5.4) and using (4.5), we find the rates for e + e annihilation in the 
form 



R E n — $2,n + °4 (Rn,2j + ^n,2j-l, 



(5.12) 



j>n-2 



for n < 5 and j < 3, where n is the number of resolved jets and R n j has i powers of 
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the logarithms (either k or A): 



R22 

R21 

-R24 
-R23 

-R26 

-R25 



-2C f kX 

\c F x 

2C 2 f k 2 X 2 



[-&o(2k + A) -3C f X]C f kX 
--CWX 3 



[2b (2K + A) + 3C f X]C 2 f k 2 X 2 



(5.13) 



R32 
R31 

R34 — 
R33 = 
R36 = 
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2C f k\ 

— -C F X 
2 

-4C F - -C A 



C F K 2 \ 2 



b (2n + X)+[^C A + 6C F -^n f )X 



C f k\ 



4C 2 F + C A C F + l -C 2 A 



1 



3 \3 



C F K X 



--C A -4C F )b (2K + \) + 
49 



-gCl - ^CaC f - 9C 2 F + ^C A n f + \c F n^j A 



C f k 2 \ 2 (5.14) 



-R44 
-R43 
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2Ci? + -Ca 



2 \2 



C F K X 



5 1 

--Ca - 3CV + g^/ 



5 
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C A + 2C F b (2K + A) + 



37 49 1 13 

^Cl + - C A C F + 9C 2 - -C A n f - -C F n f ) A 



C f k 2 X 2 (5.15) 



-R56 
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-C| + C A C F + -C 2 A 



3 \3 



C F K X' 



71 49 17 
"72^ " Y2° aCf ~ ^ + 18 CAHf + 18^ 
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One check on these results is that the DLA coefficients agree with the previous 
computation in [14] for i? 44 and R^q. A second is that the 2-jet inclusive fraction 
obtained by summing the rates is 1. 



6. Average jet multiplicity 

The average jet multiplicity, i.e. the mean number of resolved jets, as a function of 
the hard process scale Q, the angular resolution R and the minimum jet energy Er, 
provides a useful overall measure of jet activity and substructure. As indicated by 
Eq. (4.3), this quantity is obtained simply from the first derivative of the relevant 
generating function. 

6.1 Next to double leading logarithms 

Writing the average jet multiplicity in e + e~ annihilation as J\f ee , from the fixed-order 
jet rates (5.13)- (5.16) we obtain to 0(a^) 

Af ee = 2 + a s \2kX - ^XJ C F 

+ «s \ \c a k\ + 2& (2k + A) - ^C A X + jU/AJ C f kX 

/ i i i 31 \ 

+ «! ( Ys°a kX + 2 bo ° A ( 2K + A ) + Y8 nfCpX ~ 72^ A ) ° fk2X2 • (0) 

Here the terms in bo originate from including the running coupling. We see that 
these terms enhance the average jet multiplicity with respect to a fixed- coupling 
calculation. 

To perform an all-orders resummation to NDLA, we repeat the analysis of [15], 
for the inclusive algorithms instead of the exclusive kt algorithm. In terms of the 
generating functions, we have 



2 », 



du 

u=l 



2Ai q (6.2) 



du 

K(E,0 = 1+1 I dz^^P gg {z)U g {zE^) (6.3) 



u=l 

where A/"„ i9 , the average quark and gluon jet multiplicities, satisfy the equations 



f 




f 


dz^ 




y 


Je r /e 


2tt 


f 




f 


dz^ 




y 


Je r /e 


2tt 



= 1 + / -jr dz^{P gg (z)Af g (zE,a 
JUr s Je r /e 

+P qg (z)[2Af g (E,a-K(E,a]} ■ (6-4) 
In Appendix A we show that (6.3) is equivalent to the following PDE in terms 



of the logarithmic variables (4.9): 

d 2 AT g _ 

dndX \" y 4 dn 



-C F «gU,--^ , (6.5) 
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6.3 Running coupling 

Taking into account the running of a>s to next-to-leading order, we have 

(2^/C A a^\) = [1 - b {2K + A)57s + O{a 2 s )]I fey/Crf&O) . (6.15) 

Thus if we drop terms of relative order a§ the solution to (6.11) is 

N g =[l + b (2K + A)a s ] h (2^/C A a S K\) , (6.16) 

which agrees with the &o _ dependent terms in (6.1). However, for large k and/or A, 
b (2n + A)«s ~ 1 and therefore we need to take into account the running of as to all 
orders. 



6.4 Numerical solution: DLA 

Treating the running of as to all orders, as in Eq. (5.3), but still neglecting the finite 
parts of the splitting functions, we have in place of (6.11) 



d 2 K 



Ma 



(6.17) 



with c g = This PDE is not straightforward to solve explicitly. Its properties 

are discussed in Appendix B. It may however be solved numerically by discretization. 
Writing 



Af q (K = ma, X = nb) = f m>n 
N g {n = ma, X = nb) = # m>n 



(6.18) 



we have 



and 



d 2 Mg 1 
rN_i 

8k dX ab 



[Sm+l,rc+l 9m+l,n 9m,n+l 9m,n\ 



(6.19) 



CgAfg 



(2k + A + /i) 4 
+ 



9m+l,n+l 



9m+l,n 



_2{m + l)a + (n + l)b + fi 2{m + l)a + nb + fi 

9m,n+l 9m,n 



2ma + (n + l)b + \x 2ma + nb + \x 



(6.20) 



Equating these expressions, one can solve iteratively for g m+1>n+ x starting from the 
boundary values go >n = g m o = 1. 
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6.5 Numerical solution: NDLA 

To include the finite parts of the splitting functions, we may write (6.5) and (6.8) 
with equivalent precision as 

^ = c (l-d °-) N ° (621) 

where c 9;fl = Cp^/bo and 

3 , 11 rif , . 

d '=4' d ° = T2 + W>- (6 - 22) 

The PDEs (6.21) can be solved numerically by a simple extension of the method 
outlined above. For the discretized K-derivative, we use 

dU g 1 

q ~ ^ [9m+l,n+l + 9m+l,n ~ 9m,n+l ~ 9m,n\ ■ (6.23) 

We can then write the right-hand side of (6.21) as 



C 9 



,1 ~~ $g)9m+l,n+l (1 ~ &g)9m+l,n 



2(m + l)a + (n + l)b + fi 2(m + l)a + nb + ft 
(1 + 5 g )g m>n+1 (1 + 5g)g m:n 



(6.24) 



2, 2/11 n f 

^ = :^ = :kC - ( 6 - 25 ) 



2ma + (n + 1)6 + /i 2ma + nb + fi 

where 

rt. = 

a 9 a V 12 6iV 3 

and equate this to (6.19). 

Similarly, to obtain the quark jet multiplicity we write 

2 J\f q 1 

Oti OX ~ ~ob \-f rn ~^^-> n ~^^ fm+l,n fm,n+l fm,n\ j (6.26) 

equate this to (6.24) with c g , S g replaced by 

= = £ (6.27) 

to obtain the discrete equivalent of (6.5), and solve iteratively for f m +i tn +i starting 
from the boundary values /o n = / m o = 1- 



7. Monte Carlo results: e + e 

In Fig. 1 results for the average jet multiplicity from the SHERPA Monte Carlo are 
compared with the fixed-coupling DLA results (for fixed as = 0.118). The average 
jet multiplicity is shown for E cm = 200 TeV and Er = 10 GeV as a function of the jet 
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Figure 1: Average SHERPA jet multiplicity in e + e~ annihilation at 200 TeV, compared 
to the DLA formula (6.12) (solid red), and the asymptotic limit (6.14) (red dashed). Black 
and blue show SHERPA hadron and parton levels, respectively. The solid, dashed and 
dot-dashed curves show anti-fct, C/A and kt algorithms, respectively. 



radius parameter R. For these extreme values we expect the leading logarithms to 
dominate. At such a high energy, the asymptotic approximation (6.14) (red dashed) 
agrees well with the exact DLA formula (6.12) (solid red), even at R ~ 1. The Monte 
Carlo results follow the formula up to ln(l/i?) ~ 4, after which the parton-level Monte 
Carlo result saturates at the parton multiplicity (about 100 at the shower cutoff 
Qo ~ 1 GeV). Surprisingly, the hadron-level result follows the DLA perturbative 
prediction further, up to around 200 jets at ln(l/i?) ~ 6, before breaking away to 
saturate at the hadron multiplicity. We note that the quantitative agreement between 
Monte Carlo simulation, implementing the running coupling, and the analytic DLA 
result depends on the choice for the fixed coupling in the latter. For the extreme 
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cm. energy we consider, a somewhat smaller value for as can improve the agreement 
with the simulation for larger values of \n(l/R). At the same time for lower values 
of ln(l/i2), where we do not expect the DLA to fully apply, the agreement would get 
worse. 

Figure 2 shows SHERPA results on the jet multiplicity as a function of cm. 
energy, compared with the NDLA predictions for running as, derived numerically as 
explained in the previous Section. A small QCD scale A ~ 50 MeV in the NDLA 
formulae gives reasonable agreement with the Monte Carlo results, in which Ajy$ = 
180 MeV. Such a scale change is a next-to-NDLA effect and therefore within the 
uncertainties. At the smallest values of R, the Monte Carlo values are approaching 
the NDLA predictions from below, indicating further subleading effects. As in Fig. 1, 
the agreement at very small R is better at hadron level than at parton level. 
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Figure 2: Average SHERPA jet multiplicity in e + e~ annihilation as a function of centre- 
of-mass energy, compared to the NDLA result from Eqs. (6.5) and (6.8), shown in red. 
Black and blue show SHERPA hadron and parton levels, respectively. The solid, dashed 
and dot-dashed curves show anti-fcj, C/A and kt algorithms, respectively. 

It can be see from both figures that the Monte Carlo results for all three in- 
clusive algorithms, anti-fct, C/A and kt, are very similar, consistent with the lack 
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of dependence of our analytical predictions on the power p in Eqs. (2.1) and (2.2). 
The jet multiplicity is sytematically slightly lower for the k t algorithm, which we 
conjecture is due to its tendency to gather more low-momentum particles into jets 
from angles somewhat larger than the canonical jet radius R. This leads to a smaller 
number of jets for a given final-state multiplicity. This higher susceptibility of kt jets 
to radiation was discussed in [16] where also a first-order estimate for the resulting 
effective area of k t and C/A jets was derived. Qualtitatively we observe the same 
behaviour as for the passive jet area discussed there: the effect increases for higher 
jet energies and larger nominal R. 

8. Monte Carlo results: pp 

The motivation for studying the inclusive family of jet algorithms is the similarity to 
phenomenologically relevant hadron collider algorithms. For jet rates, and therefore 
also the mean number of jets, the factorization of the PDF in the initial state holds 
at the double leading logarithmic order. In this section we compare the predictions of 
our generating functions with Monte Carlo calculations and determine the reliability 
of our resummation in the LHC context. 



Mean jet multiplicity at 8 TeV LHC (Sherpa MC) 




H T /2 

Figure 3: The proton-proton average jet multiplicity in inclusive di-jet production as a 
function of Ht/2. Compared is the DLA analytic prediction, shown in red, with parton 
shower results obtained with SHERPA. The solid, dashed and dot-dashed curves show 
anti-/c(, C/A and kt algorithms, respectively. 
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Figure 3 shows results for pp collisions at the current LHC energy of 8 TeV. Here 
the average multiplicity of jets with p? > En = 20 GeV is plotted as a function of 
Ht/2, where 



the scalar sum of all jet transverse momenta. This choice corresponds closely with 
the initial parton transverse momenta in the QCD 2 — > 2 hard scattering. Each 
point in the analytic result is weighted by the partonic fraction of final-state quarks 
versus gluons. In other words, the analytic result is 



where c q + c g = 1 and are determined from the proportions of final-state jets in the 
contributing 2—7-2 hard subprocesses. The considered values for R are the same as 
in Fig. 2. The multiplicity levels off and even decreases at high values of the primary 
parton p?, owing to the transition from gluon to quark jets. Overall we observe a 
good agreement of the analytic estimate with the SHERPA simulation. As before 
we use a fixed coupling of as = 0.118 in the DLA calculation. In particular for large 
values of Ht/2 a reduced value could be more appropriate, improving the agreement 
with simulation. Notably, the dependence on the jet algorithm for the simulated 
results is rather mild. As observed for e + e~ collisions in the above, the k t algorithm 
tends to produce slightly less jets for a given final-state multiplicity. 

9. Jet scaling patterns 

The work in [17-19] discussed the potential for extrapolating jet rates based on 
universal scaling patterns. These patterns are most easily classified in terms of 
the ratio of exclusive jet rates R( n +i)/n = cr n+ \/a n . In the Durham (exclusive k t ) 
algorithm, it was found that for low multiplicity n < (nj ets ), emissions are essentially 
Poisson-like so that Rr n+ iy n ~ (n + 1) _1 . The tail of the multiplicity distribution 
then produces dominantly staircase or geometric scaling where Rt n+ iy n ~ constant. 
This regime is driven by the fractal nature of QCD radiation in the gluon dominated 
limit. 

We know from previous work that the expected scaling patterns of jets can de- 
pend dramatically on the jet algorithm. One example of this is the JADE algorithm, 
where the non-exponentiation of the primary emissions precludes the Poisson extrap- 
olation even in the pseudo-abelian limit [20]. In this section we would like to address 
scaling in the inclusive generalized k t class of algorithms. This extends the results 
in [17] and strengthens the case for investigation at hadron colliders. 




(8.1) 



(rijj) = 2(c q (n q ) +c g (n g )), 



(8.2) 
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Figure 4: Jet ratios from SHERPA compared with the Poisson extrapolation from the 
first bin. For smaller jet sizes the Poisson behaviour breaks down as discussed in the text. 
Jet sizes here correspond to approximate R values of 0.61 (left) and 0.22 (right). 



9.1 Poisson breaking components 

With the leading logarithmic coefficients from Eqs. (5. 13)- (5. 16) it is easy to make 
some first statements about scaling in the generalized algorithm. It is clear from the 
structure of the coefficients that in the limit Ca — > a perfect Poisson distribution 
emerges. Now a simple comparison between the generalized and Durham algorithms 
is the relative size of the Poisson breaking components in the lower multiplicity 
rates, for example the 2-gluon correlated emission contribution to the 4-jet rate 
(5.15). For the double-leading logarithmic coefficients to the 4-jet rate, -R44, we 
find Q D -ham _ 2C 2 p + (i/^CaCf and Cfr ralizcd „ 2C| + (1/2)C A C F using the 
normalization of this work. The lowest-order Poisson breaking term is relatively 
larger in the generalized class of algorithms. We would thus expect the onset of 
staircase (geometric) scaling to come about even for smaller values of the logarithm, 
and to better match the staircase behaviour at lower multiplicity. Evidence of this 
may be found in the Z + jets analysis presented in [17]. We present in table 1 the 
relative sizes of the Poisson breaking terms in the DLA expanding coefficients for 
the two algorithms compared with the idealized Poisson and staircase predictions. 
In this table we compute the 6-jet rate in the generalized algorithm by expanding 
the resummed jet rates from Eq. (5.11) to 0(af) which entirely determines the DLA 
resolved component. 
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Poisson 


Generalized k t 


Durham 


Staircase 


-R4/3/ Rz/2 


0.50 


0.781 


0.688 


1 


-R5/4/ -R4/3 


0.67 


0.906 


0.868 


1 


^6/5/ -^5/4 


0.75 


0.923 


0.932 


1 



Table 1: Ratio of successive ratios for the generalized kt and Durham jet rates compared 
with the idealized Poisson and staircase expectations. 

9.2 independence of scaling 

A second question we wish to answer in this section is how the idealized scaling 
patterns depend on the jet radius parameter R in the generalized kt algorithm. The 
additional handle provided by the separation of the angular and energy regulation 
allows us to probe an effect unbeknownst in the Durham algorithm. At the double 
leading logarithmic level, the resummed rates in the general algorithm are invariant 
under exchange of «-<->■ A. Decreasing the size of the jet merely increases the overall 
logarithm. 

In the simulation however, although we increase the overall logarithm when we 
require smaller resolved jets, we also change the relative contributions of the primary 
and secondary contributions due to kinematics. This effect is present even at the 
level of the ordered 2-gluon emission matrix element [21]. 

Using the parton shower we find that larger jet sizes dramatically increases the 
goodness of the fit with respect to the Poisson hypothesis. This has a clear inter- 
pretation. Correlated emissions in these events are predominantly intra-jet evolution 
for large jet radius. In addition, this effect is larger than that due to the breaking of 
K <H- A symmetry at NDLA. 

To analyze this further we investigate the observable #34, defined as 

a _ P3 • P4 /„ 1 s 

cos6»34 = 1 — r, (9.1) 

|Ps||P4| 

in reconstructed 4-jet events. The third and fourth hardest jets overwhelmingly 
originate from emitted gluons. As we suspect that in the simulation secondary Pois- 
son breaking emissions are on average closer in # 34 , this observable gives another 
estimate on the size of these breaking effects. Comparing the prediction of pseudo- 
abelian QCD (in practice the parton shower with the g — >■ gg and g — >• qq branching 
off) we see in Fig. 5 that the jet size has a large effect on the relative size of the 
correlated emission component. In fact it appears that roughly R ~ 0.4 — 0.5 is 
the smallest value of the radius where the primary emissions are still the dominant 
contribution to the 4-jet rate. 
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Angle between 3rd and 4th hardest jet (Sherpa MC) 




0.5 1 1.5 2 2.5 3 

#34 



Figure 5: Differential cross-section with respect to the angle #34 in exclusive 4-jet events. 
The solid curves in both cases represents the parton shower while the dashed curves have 
the gluon splitting function turned off (both to quarks and gluons). The large contribution 
of subsequent splittings between the red curves (small jet sizes) prevents a valid Poisson 
extrapolation as seen in Fig. 4. 

10. Conclusions 

In the present paper we have derived the generating functions of jet rates for the 
inclusive generalized k t jet algorithms, valid in the next-to-double log approximation 
(NDLA), and used them to compute jet rates and average jet multiplicities as func- 
tions of the jet radius parameter R and the minimum jet energy Er. At this level of 
precision, the results are independent of the power p that distinguishes the inclusive 
k t , Cambridge/ Aachen and anti-fcf algorithms. 

The analytical results on e + e~ annihilation are in broad agreement with those 
from the SHERPA Monte Carlo event generator, including the weak p-dependence. 
Surprisingly, the hadron-level Monte Carlo results follow the analytical predictions 
down to very small jet radii, well beyond the range of the perturbative parton shower, 
indicating that the cluster hadronization model in SHERPA is smoothly matched 
to perturbation theory. Analytical predictions for pp collisions at the current LHC 
energy also agree fairly well with SHERPA Monte Carlo results. 

The SHERPA Monte Carlo results and generating function in the generalized 
k t algorithm were also used to study the transition from Poisson-like to staircase-like 
behaviour in the (n+ l)/n jet ratios. The relatively larger non-abelian terms in the 
DLA e + e~ jet rates, compared to the Durham algorithm, indicate a transition at 
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lower values of n. Furthermore, the jet radius was found to have a significant impact 
on scaling, where smaller jet sizes receive larger secondary contributions. This feature 
may be relevant for jet activity studies in jet substructure analyses. 

These results represent the most detailed comparison to date between analytical 
and Monte Carlo predictions for the inclusive jet algorithms in current use at the 
LHC. Information about the dependence of jet rates and multiplicities at high jet 
energies and relatively small radii are particularly important for studies of boosted 
jets and jet substructure, which play an increasingly important role in searches for 
new physics at the LHC. Our leading-logarithmic predictions are relevant to such 
studies, and their fairly good agreement with Monte Carlo results is encouraging. 
Possible improvements such as optimized scale choices and empirical subleading con- 
tributions would be worth exploring. Studies of other features of the jet multiplicity 
distribution along these lines could also be useful. 
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A. Derivation of the PDE for the average jet multiplicity 

Differentiating Eq. (6.3) with respect to £ and using (3.3), we have to NDLA 



But at the same level of precision we have (remembering the £ dependence of as) 



d Mq-l 
<9£ «s 



C F 




boas-\n.z)M g {zE,£) 




(A.l) 



d M q -l _ 1 dM q 



+ h{M q - 1) 



(A.2) 



<9£ a s a s <9£ 



and 



C F a s 




■M g {zE,£')=M q -\, 



(A.3) 



so that 




(A.4) 
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Differentiating with respect to E and using (3.3) again 
d i dK 



= Cf ^ + J Cf I dz WW,® -Af 9 (E,0\ 
a s at, 4 Jer/e 

+2C F b a s / —N g {zE,£) . (A.5) 

J En IE Z 



Now to the required precision 



f 1 dAf r 1 

/ dz[M g {zE^)-M g {E^)] = E—^ dz\nz = -E 
Je r ie vE Jo 

and 

dz if/ t i >, \ M, 



dE 



(A.6) 



'e r /e 

while (remembering the E dependence of as 



f dz dU 
2C F b a s / -Af g (zE, = 2b ^ (A.7) 



oE as at, as oE o£, ot, 



Thus we get the PDE 

^ d dJ\f q „ ( \ r 3 rn dAf g \ 

E dE^ = c ^{^-i E ^i) • (A - 9) 

which gives (6.5) in terms of the logarithmic variables (4.9). 

B. Properties of the PDE for the average jet multiplicity 

The PDE (6.17) is not so straightforward to solve explicitly. If we consider changing 
to variables corresponding to the sum and difference of k and A, we find that the PDE 
is separable (one obtains ODEs that are equivalent to those describing the classical 
mechanics of the harmonic oscillator and the quantum mechanics of the Coulomb 
potential), but that the boundary conditions are not. As a result, one may prove that 
one cannot express the boundary conditions as a Fourier-Bessel series of orthogonal 
functions of Sturm-Liouville type. 

On reflection, seeking an explicit solution in this way is perhaps ambitious, given 
that there is no guarantee that solutions to such a PDE can be written in terms of 
special functions. We thus proceed to an analysis of a rather different nature. This 
analysis will enable us to (i) establish that a solution exists; (ii) give a variety of 
infinite montonic series of upper and lower bounds on that solution; and (iii) provide 
explicit series solutions that enable, e.g., the asymptotic behaviour of the solution 
to be found in simple, closed form. 
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To proceed, it is useful to define x = ^-(2k + |), y — y(A+ |), in terms of which 

M xy = (B.l) 

x + y 

We may recast the PDE in the integral forms 

Z'C^ < b - 2 » 

where a = ^p, or 



dy J a x' + y'' 



v 



" = 1 + J.«J.'"*Tf (B ' 3) 

The first of these forms, by the way, makes clear that a physically-acceptable solution, 
Af, if it exists, is a monotonically increasing function of x, V y, and of y, V x. (Proof: 
The mean number of jets must be > 0. Thus Af y > 0; Af x > follows by symmetry 
inifi y.) 

The basic idea behind our analysis will be to first identify related PDEs for which 
we can find explicit solutions whose properties (such as their existence, continuity, 
&c.) can be checked 'by hand'. We then use these solutions as a crutch to derive 
properties (for example, the existence) of the solution of the original PDE. 

Our first lemma is as follows. Suppose M°(x,y) is a solution of the PDE 

M° xy = a(x,y)M° (B.4) 

subject to the same boundary conditions as Af. Suppose furthermore that a(x,y) > 
■^y- almost everywhere 1 in the domain D = {(x,y)\x,y > a}. Then 

rv r x /wo 

M\x,y) = 1+1 dy' dJ—— < M°, V (x,y) e D. (B.5) 

J a J a X+y 

( For brevity, we denote a double integral of the type that appears on the RHS by 
J henceforth.) Proof of this follows immediately from the fact that 1 + J < 
1 + J aA4° = A4° everywhere in D. Continuing in this vein, we define an infinite 
sequence whose elements are 

/M i 
^- . (B.6) 
x + y 

Now 

M + -M = / (B.7) 



x + y 

so M i - M 1 - 1 < V(x,y) e D M l+l - M l < V(x,y) e D and proof that 

this sequence is montonically decreasing follows by induction, since we have already 
proven that A4 1 < A4° everywhere in D. 



1 We use the usual language of measure theory. 
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We next prove that the sequence {.M* } is bounded below by zero, if Ai° > 
(which we can verify by hand given an explicit Ai°), and hence converges. To do 
so, we simply note that positivity of Ai° implies positivity of M. % V i by induction, 
given (B.6) and M° > 0. 

We now prove that the sequence {M 1 } converges to a solution, Jf of the original 
PDE. Indeed 

Jf = lim M = lim M i+1 = 1 + lim / = 1 + / lim = 1 + 



x + y J i->-oo x + y J x + y 

(B.8) 

where, in the penultimate step we used the Theorem of Monotone Convergence of 
Lebesgue integration. 2 

We shall call Jf a supersolution of (B.l). Evidently, given any suitable a and Ai°, 
we can use the above results to prove that a solution exists and to find an infinite 
series of monotonically decreasing upper bounds on it. This raises the question 
of uniqueness, however: different starting points, (a,Ai°), may lead to different 
supersolutions. 

It is straightforward to prove the following weak version of uniqueness, for se- 
quences that are 'nested' in the following sense. Suppose two sequences and 
{Ai' 1 } converge on supersolutions Jf and AT, respectively. If 3 i,j,k,l such that 
elements M. 1 > Ai'i and Ad < At' 1 almost everywhere, then Jf = Jf . Proof: 
M i > M' j M* > Jf. But M i - Jf > f = M i+1 -Jf > 

— — — J x+y — 

Jf > Jf'. Similarly, M k < M' 1 =^ Jf < Jf =^ Jf = Jf , QED. Unfortunately, 
this weak version of uniqueness does not even imply that any two supersolutions 
coincide, since the sequences that lead to them may be wholly contained in distinct 
intervals, or, even if they are nested, an element of one sequence need not exceed an 
element of the other almost everywhere. But this weak version is still useful where it 
applies, in that it may give us a more rapidly converging sequence of upper bounds 
on a solution. 

Let us now discuss lower bounds. Similar arguments to those just given show 
that if we instead started with a solution, C°(x,y) to a PDE of the form 

C° xy = P(x,y)C° } (B.9) 

with f3 < -j- almost everywhere on D, then we shall obtain an infinite, monotonically 

x-\-y 

increasing sequence. We can, furthermore, prove that this sequence is bounded above 
and hence converges, provided we can show by hand that CP < A4°, for some Ai° 



2 This Theorem requires that the M 1 be measurable, but this follows from the fact M°, and 
hence M 1 are C°. Note, however, that we have not proven that Af itself is C°, let alone C 1 or C 2 . 
Strictly speaking, therefore, we have proven that J\f is a solution of the integral equation, which 
is what we started with, rather than the PDE. We shall in any case ignore this subtlety in what 
follows. 
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that converges to a supersolution A/". To wit, suppose we have shown that £° < 
It then follows by induction that D < M\ since C t+1 — M' l+l = J Cl ~^' , and 
thus that linij^oo C 1 < A/*. Again, by the Theorem of Monotone Convergence, {£*} 
converges to a solution of (B.l) that we call a subsolution and denote by M_. By a 
straightforward generalization of the result just proven, a subsolution is less than or 
equal to any supersolution for which it can be shown that any one element in the 
sequence defining the former is less than or equal to (almost everywhere) any element 
in the sequence defining the latter. 

Yet again, this result is insufficient to establish uniqueness, but it does provide 
a way to obtain both upper and lower bounds on a given solution. Obviously, if we 
can prove (or assume) uniqueness, then all of the aforementioned super- and sub- 
solutions coincide, meaning that any of the aforementioned sequences can be used to 
bound the solution to arbitrary, known precision. 

We can prove uniqueness of the solution heuristically out to any finite (x, y) by 
explicit construction. Given J\f — 1 on the boundary, we divide the intervals [a, x] 
and [a, y] into regions of size Sx and Sy and find the unique solution Af(a + Sx,y) ~ 
l+5x -£^7, with a similar result for A"(x, a+5y). We then extrapolate to the region 
[a + 5x, a + 25x] and so on. Given that we have proven the existence of the solution, 
we do not need to worry that the extrapolations in the x and y directions might 
not coincide. To make this proof rigorous outside a neighbourhood of the boundary, 
we would then have to take the limit 8x, 5y — > 0. We will content ourselves with 
assuming that the solution is unique. 

As we stressed above, all of our results on super- and sub-solutions are contingent 
on showing the existence of solutions to the equations (B.4) and (B.9), for suitable a 
and (3, satisfying the BCs, and showing that they have the desired properties. We do 
this by supplying explicit solutions for functions a, (3 of the combination (x— a)(y— a), 
for which the PDE reduces to an ODE, and for which the boundary conditions reduce 
to a single boundary condition at {x — a)(y — a) = 0. 

To obtain lower bounds, once may start, for example, with the solution C° = 1 
of the PDE C° xy = 0, since > everywhere in D. We thus obtain 

C 1 > 1 + (x + y) ln(x + y) - (x + a) ln(z + a) - (a + y) In (a + y) + 2a In 2a, (B.10) 

with subsequent D given in terms of polylogarithms. At fixed y and large x, for 
example, we deduce that C 1 ~ (y — a) In a;. More generally, C n ~ ^^~Stw ^ =>• 



To find an upper bound, we can use a == ^ which evidently exceeds 
throughout D. This yields an initial upper bound of the form 



o(V - a)\nx) 



y/2n(y-a)logx 



M 







I Q {y/2a{x - a){y - a)). 



(B.ll) 
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A more stringent upper bound may be obtained from a = — , 1 > -j-, with 

2y/ (x-a)(y-a) X+ V 

solution 

M° = I (2y/2(x-a)*(y-a)*). (B.12) 
For a better lower bound, we can, provided a > 1, solve 

2 

£ ^ = (^F £ °' (R13) 

whose solution is 

£° = 7 + °' (B.14) 

It is soothing to verify explicitly that the resulting subsolution has the same asymp- 
totic behaviour, at large x and fixed y, as the subsolution starting from L° — 1 
derived above. Indeed, in this limit, C° ~ -, C 1 ~ y ~ a logx, and so on, with 

C n _ ( y -a)"(i/+na) io^x At k we thug find "£n _ {y^a£io^ exactly as before. 

o(n+l)! n\ ° ' n! n! ' J 
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